###### Maintext: Table 2
###### Treatment Effects on Punishment, Ownership, and Ladder Position
gc(); rm(list = ls()); set.seed(12345)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) # Note: if you are not using R Studio this command will not work, set WD to source file location manually

############### SCRIPT SUMMARY #####################

packages <- c("dplyr", "ggplot2", "texreg")
lapply(packages, require, character.only = T)
source("functions.R")

### Load Datasets

dfMerge <- read.csv("data/cleaned/uganda_jun17_cleaned.csv", header=T, stringsAsFactors =F)
dfMerge <- dfMerge[!dfMerge$game_id_exp %in% c("Basic Aid", "Aid Ownership"), ]

########################################## MISC PREP AND VARGEN ######################################

## Re-leveling treatment variable 
print(table(dfMerge$treatment_pooled))
newLabs <- c("Hidden VAT", "Direct Tax", "Visible VAT", "Windfall") 
origLabs <- c("Strong Hidden VAT", "Direct Tax", "Weak Hidden VAT", "Windfall")
treatCor <- setNames(newLabs, origLabs)
dfMerge$treatment <- treatCor[match(dfMerge$treatment_pooled, 
                                           names(treatCor))] 
table(dfMerge$treatment)
dfExp <- dfMerge
dfExp <- dfExp[dfExp$role == "C" & dfExp$round != 1 & 
                 !is.na(dfExp$treatment), ] # First round is a practice round

## Determine which columns to take from the model frames

colsInc <- c("outcome", "comparison", "estimate", "std.error", 
             "conf.low", "conf.high", "p.value", "n")

######################## PANEL A: THRESHOLD AS DV

require(estimatr)

refCat <- c("Visible VAT", "Hidden VAT", "Direct Tax", "Windfall")
dfComp <- c()

for(i in 1:length(refCat)){
  
  depVar <- c("sbj_thresh")
  dfUsing <- dfExp
  dfUsing$treatment <- relevel(factor(dfUsing$treatment), 
                                      ref = refCat[i])
  expVars <- c("treatment", "as.factor(enum_id)", 
               "item", "as.factor(round)")
  controlVars <- c("age", "edu", "poverty_avg", "services_avg", "female")
  
  modEq <- as.formula(make_equation(depVar, c(expVars, controlVars)))
  modFit <- lm_robust(modEq, data = dfUsing, 
                          clusters = sbj_unique,
                          se_type = "CR2")
  
  dfMod <- tidy(modFit)
  dfMod$reference <- refCat[i]
  dfMod$term <- gsub("treatment", "", dfMod$term)
  dfMod$n <- modFit$nobs
  dfMod <- dfMod %>% filter(term %in% newLabs)
  dfComp <- rbind(dfComp, dfMod)
}

## Subset to only the comparisons you want

dfComp$comparison <- paste(dfComp$term, "-", dfComp$reference)
compInc <- c("Hidden VAT - Visible VAT", "Visible VAT - Direct Tax", 
             "Hidden VAT - Windfall", "Visible VAT - Hidden VAT", 
             "Direct Tax - Hidden VAT")
(dfThresh <- dfComp[dfComp$comparison %in% compInc, 
                 colsInc])


#################### PANEL B: OWNERSHIP AS DV 

dfUsing <- dfExp[dfExp$round == 2, ] # Subset to single round since it's ownership

require(estimatr)

refCat <- c("Visible VAT", "Hidden VAT", "Direct Tax", "Windfall")
dfComp <- c()

for(i in 1:length(refCat)){
  depVar <- c("strong_ownership")
  dfUsing$treatment <- relevel(factor(dfUsing$treatment), 
                               ref = refCat[i])
  expVars <- c("treatment", "as.factor(enum_id)")
  controlVars <- c("age", "edu", "poverty_avg", "services_avg", "female")
  
  modEq <- as.formula(make_equation(depVar, c(expVars, controlVars)))
  modFit <- lm_robust(modEq, data = dfUsing, 
                      se_type = "HC3")
  
  dfMod <- tidy(modFit)
  dfMod$reference <- refCat[i]
  dfMod$term <- gsub("treatment", "", dfMod$term)
  dfMod$n <- modFit$nobs
  dfMod <- dfMod %>% filter(term %in% newLabs)
  dfComp <- rbind(dfComp, dfMod)
}

## Subset to only the comparisons you want

dfComp$comparison <- paste(dfComp$term, "-", dfComp$reference)
(dfOwn <- dfComp[dfComp$comparison %in% compInc, 
                 colsInc])


######################## PANEL C: LADDER AS DV

require(estimatr)

refCat <- c("Visible VAT", "Hidden VAT", "Direct Tax", "Windfall")
dfComp <- c()

for(i in 1:length(refCat)){
  
  depVar <- c("l1")
  dfUsing <- dfExp
  dfUsing$treatment <- relevel(factor(dfUsing$treatment), 
                               ref = refCat[i])
  expVars <- c("treatment", "as.factor(enum_id)", 
               "item", "as.factor(round)")
  controlVars <- c("age", "edu", "poverty_avg", "services_avg", "female")
  
  modEq <- as.formula(make_equation(depVar, c(expVars, controlVars)))
  modFit <- lm_robust(modEq, data = dfUsing, 
                      clusters = sbj_unique,
                      se_type = "CR2")
  
  dfMod <- tidy(modFit)
  dfMod$reference <- refCat[i]
  dfMod$term <- gsub("treatment", "", dfMod$term)
  dfMod$n <- modFit$nobs
  dfMod <- dfMod %>% filter(term %in% newLabs)
  dfComp <- rbind(dfComp, dfMod)
}

## Subset to only the comparisons you want

dfComp$comparison <- paste(dfComp$term, "-", dfComp$reference)
(dfLoss <- dfComp[dfComp$comparison %in% compInc, 
                    colsInc])


### Combine into single frame:

dfAll <- as.data.frame(rbind(dfThresh, dfOwn, dfLoss))
roundVars <- c("estimate", "conf.low", "conf.high", "p.value")
dfAll[, roundVars] <- round(dfAll[, roundVars], digits = 3)

## Subset to only the comparisons you want for each

useComps <- c("Visible VAT - Hidden VAT", 
              "Visible VAT - Direct Tax", 
              "Hidden VAT - Windfall")
useCols <- c("outcome", "comparison", "estimate", 
             "std.error", "p.value", "n")
(dfUse <- dfAll[dfAll$comparison %in% useComps, useCols])

texObj <- prepTexreg(dfUse, names.col = "comparison")
texreg(texObj$obj, stars = c(0.01, 0.05, 0.1), 
          custom.model.names = c("Threshold", "Ownership", "Loss"))

save.image(file = "results/tab_2.RData")


